set.seed(1234)
library(splines)
library(ggplot2)
library(viridis)
library(patchwork)
library(hexbin)
library(reshape2)
library(ResistPhy)
library(ape)
library(posterior)
library(cmdstanr)
library(bayesplot)
library(latex2exp)
run_mcmc <- F
We Restrict ourselves to a 20 year span Load usage data and simulation script
usage_df <- simulated_usage()
source("simu_phylo.R")
q_u <- 1.2
q_t <- -2.7
n_tip <- 200
gamma_sus <- 1/60.0*365.0
gamma_res_u <- gamma_sus + q_u
gamma_res_t <- gamma_sus + q_t
t0 <- min(usage_df$time)
tmax <- max(usage_df$time)
sim <- sample_phylo(gamma_sus, gamma_res_u, gamma_res_t, n_tip, 123456, 1234567)
df <- sim$traj
phys <- sim$phys
mr <- sim$most_recent_samp
N_0 <- 3e6
incidence <- 1.5e5 / N_0
mu <- incidence+gamma_sus
time_grid <- seq(from=t0, to=tmax, by=1e-4)
t_change <- (t0 + floor(0.33*(tmax-t0)))
beta_func <- function (t) if (t < t_change) 1.02*mu else mu*(0.98 + 0.06*(t-t_change)/(tmax-t_change))
N_func <- function(t) N_0*exp(log(2)*(t-t0)/(tmax-t0))
usg_func <- yearly_usg_stepfunc(usage_df$usage, usage_df$time)
beta_vals <- sapply(time_grid, beta_func)
beta_df <- data.frame(time=time_grid, beta=beta_vals)
usg_vals <- sapply(time_grid, usg_func)
usg_df <- data.frame(time=time_grid, usg=usg_vals)
N_vals <- sapply(time_grid, N_func)
N_df <- data.frame(time=time_grid, N=N_vals)
N_plt <- ggplot(N_df, aes(x=time, y=N))+
geom_line() +
xlim(t0,tmax) +
labs(y=TeX(paste0("$","N(t)","$")))+
theme_minimal()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
beta_plt <- ggplot(beta_df, aes(x=time, y=beta))+
geom_line()+
xlim(t0,tmax) +
labs(y=TeX(paste0("$","\\beta(t)","$")))+
theme_minimal()+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
usg_plt <- ggplot(usg_df, aes(x=time, y=usg))+
geom_line()+
xlim(t0,tmax) +
labs(x="Time",y=TeX(paste0("$","u(t)","$")))+
theme_minimal()
f_plt <- N_plt/beta_plt/usg_plt
pdf("Figures/functions.pdf",4,6)
f_plt
dev.off()
#> png
#> 2
f_plt
usg_fun <- yearly_usg_stepfunc(usage_df$usage, usage_df$time)
ggplot(data.frame(time=time_grid, usage=sapply(time_grid, usg_fun)), aes(x=time, y=usage))+geom_line()+theme_minimal()
plot(ggplot(df, aes(x=t))+geom_line(aes(y=I_2),color="blue") +geom_line(aes(y=I_3),color="red"))
Preview the phylogenies run BNPR to see how that looks
plot(phys[[1]])
axisPhylo()
plot(phys[[2]])
axisPhylo()
Run MCMC
out <- infer_costs2(phys,
mr,
usage_df$usage,
usage_df$time,
t0,
tmax,
1/60.0,
365,
n_iter=2000,
n_warmup=2000,
model="nodecay",
K=60,
L=6.5,
gamma_log_sd=0.3,
stan_control=list(adapt_delta=.99,
max_treedepth=13,
parallel_chains=4,
chains=4,
refresh=10
))
saveRDS(out, "~/simu_rbf.rds")
out <- readRDS("~/simu_rbf.rds")
print(out$converged)
#> [1] TRUE
source("Figures/figure1.R")
pdf("Figures/figure1.pdf",6,6)
plot(fig1)
dev.off()
#> png
#> 2
fig1
source("Figures/figure2.R")
#> Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
#> ℹ Please use the `linewidth` argument instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
#> Warning: `aes_()` was deprecated in ggplot2 3.0.0.
#> ℹ Please use tidy evaluation idioms with `aes()`
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
pdf("Figures/figure2.pdf",6,6)
plot(fig2)
#> Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
#> ℹ Please use `after_stat(density)` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
dev.off()
#> png
#> 2
fig2
s1 <- plot_traces(out)
pdf("Figures/s1.pdf",10,10)
plot(s1)
dev.off()
#> png
#> 2
s1
s2a <- plot_ppcheck_At(out,1)
s2b <- plot_ppcheck_At(out,2)
pdf("Figures/s2.pdf",6,6)
plot(s2a)
plot(s2b)
dev.off()
#> png
#> 2
s2a
s2b
s3 <- plot_hyperpar_pairs(out)
pdf("Figures/s3.pdf",6,6)
plot(s3)
dev.off()
#> png
#> 2
s3